Groupe:
1 - Onorato, Claudia, Matricule: 1845448
2 - Harvey, William, Matricule: 1851388
Ce laboratoire a pour objectif de manipuler et analyser la manipulation et le filtrage d'images dans le domaine temporel et fréquentiel.
Remise:
La date de remise est le lundi 2 Novembre à 23h55. Une pénalité de 3 points par jour sera appliquée lors d'un retard.
Documents à remettre :
Les exercices doivent être codés dans ce fichier TP.ipynb. Les réponses aux questions doivent être incluses dans le code sous forme de commentaires ou dans des cellules dédiées (Markdown ou text). Les exercices doivent être séparés par des cellules, suivant le template fourni. Vous devez bien identifier chaque exercice et sous-question, et bien commenter le code. Veuillez nommer vos variables de manière explicite et assurez-vous que toutes les figures soient lisibles.
Créer un fichier de rendu html (Fichier -> Télécharger au format... -> HTML) de votre code et de vos graphiques. Veuillez remettre tous vos fichiers (.ipynb, html et autres) dans un seul fichier zip et nommez ce fichier selon vos matricules (Mat1_Mat2.zip).
Une pénalité de 3 points sera appliquée si ces consignes ne sont pas respectées.
Première séance: Filtrage spatial d'image¶
Note: Lors de la conception de vos algorithmes, portez attention aux types de vos données (np.uint8, np.float32, etc.). Lors de la manipulation des images, il sera probablement nécessaire de passer en np.float32 pour faire vos calculs puis de revenir en np.uint8 pour afficher vos images
Note 2: Nous vous suggérons fortement d'utiliser la fonction imshow de la manière suivante:
plt.imshow(mon_img, vmin=0, vmax=255)
# ou bien
ax.imshow(mon_img, vmin=0, vmax=255)
Cela précise à matlplolib que toutes les valeurs de l'image inférieures ou égales à zéros doivent être affichées en noires et celles supérieures ou égales à 255 en blancs (ce qui l'intervalle standard d'image codée en uint8). Si vous travaillez avec des images dont les valeurs sont comprises entre 0 et 1, vous aurez à modifier l'intervalle en conséquences.
import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['image.cmap'] = 'gray' # Choix de la color map par défaut, ne pas modifier
import matplotlib
matplotlib.rcParams['figure.figsize'] = (15.0, 10.0) # Taille des figures par défaut, peut-être modifié au cas par cas.
Exercice 1 - Réhaussement d'image¶
Le but de cet exercice est de corriger les effets d'une image fortement altérée. Pour cela, vous allez implémenter un certain nombre de filtres basiques pour corriger les différents défauts.
Question 1 (1.5 points)¶
Chargez l'image TempsModernes.jpeg en utilisant la fonction:
img = cv2.imread(filename, cv2.IMREAD_GRAYSCALE)
(a) Affichez l'image. Comme vous pouvez le constater, l'image est très sombre et peu contrastée. Affichez à côté de l'image son histogramme en utilisant la fonction fournie dans la cellule ci-dessous.
(b) Nous allons corriger son exposition et son contraste en implémentant l'égalisation d'histogramme telle que vue en cours. Complétez le code suivant:
def equalize_histogram(img):
imegale = np.zeros(img.shape)
hist, bins = np.histogram(img, 256,[0,256])
hist = hist.astype(np.float32)
...
Affichez l'image égalisée et son histogramme. Que constatez-vous?
FILENAME = 'TempsModernes.jpeg'
img = cv2.imread(FILENAME, cv2.IMREAD_GRAYSCALE)
def plot_histogram(img, ax=None):
"""
param img: image au format np.ndarray, dont les valeurs sont comprises entre 0 et 255
param ax: (optionnel) axe matplotlib sur lequel tracer l'histogramme.
"""
hist, bins = np.histogram(img, 256,[0,256])
if ax is not None:
ax.bar(bins[:-1], hist)
else:
plt.bar(bins[:-1], hist)
plt.show()
def plot_image_with_histogram(img, description):
img = img.astype(np.uint8)
fig, ax = plt.subplots(ncols=2, figsize=(20, 5))
ax[0].imshow(img, vmin=0, vmax=255)
ax[0].set_title(f'{description}')
plot_histogram(img, ax=ax[1])
ax[1].set_title(f"Histogramme de {description}");
fig.show()
plot_image_with_histogram(img, description="image originale")
def equalize_histogram(img):
imegale = np.zeros(img.shape)
normalized_hist, bins = np.histogram(img, 256,[0,256], density=True)
normalized_hist = normalized_hist.astype(np.float32)
normalized_cum_hist = np.cumsum(normalized_hist)
for i in range(imegale.shape[0]):
for j in range(imegale.shape[1]):
# since img is in uint8, the values are already equivalent to the hist bins.
index = np.argwhere(bins == img[i, j])[0][0]
imegale[i, j] = 255 * normalized_cum_hist[index]
return imegale.astype(np.uint8)
def equalize_histogram_vectorized(img):
imegale = np.zeros(img.shape)
normalized_hist, bins = np.histogram(img, 256,[0,256], density=True)
normalized_hist = normalized_hist.astype(np.float32)
normalized_cum_hist = np.cumsum(normalized_hist)
# Vectorized operation; instead of iterating over each pixel of the original image, we use the facts that:
# - The values contained in the original images directly maps to the indexes of the bins, because the possible values range from 0 to 255
# - The bins values and indexes are equal
# the operation then takes a lot less time.
assert set(np.unique(img)).issubset(set(np.unique(bins))), "the values contained in the original images has to directly map to the values of the bins"
assert (bins == list(range(len(bins)))).all(), "bins values and indexes are equal"
imegale = 255 * normalized_cum_hist[img]
return imegale.astype(np.uint8)
%%time
equalized_img = equalize_histogram(img)
plot_image_with_histogram(equalized_img, description="image égalisée")
%%time
equalized_img_2 = equalize_histogram_vectorized(img)
plot_image_with_histogram(equalized_img_2, description="image égalisée")
Réponse:
En faisant égalisation l'histogramme de l'image, nous pouvons remarquer que le contraste a grandement été augmenté dans l'image. Contrairement à l'image originale, qui est très sombre, on peut maintenant voir dans l'image tant des coins sombres que des coins lumineux. De la même manière, on voit que l'histogramme, plutôt que d'avoir des valeurs majoritairement concentré entre 0 et 100, est maintenant étalé de 0 à 255.
Question 2 (1 point)¶
Quel type de bruit est semble particulièrement prédominant dans l'image égalisée? Expliquez pourquoi un filtre médian est adapté pour le corriger. Vérifiez que c'est bien le cas en filtrant l'image avec un filtre de taille 5.
Réponse: Le bruit prédominant dans l'image est un bruit poivre et sel. Ce bruit, qui est de haute fréquence, comme l'on passe d'une couleur moyenne à soit noir ou blanc, est filtré par des filtres passe-bas. Le filtre médian, qui est un filtre passe-bas, est bien adapté pour filtrer le bruit poivre et sel, comme il préserve bien les contours, contrairement au filte moyenneur par exemple.
img_median_filtered = cv2.medianBlur(equalized_img_2, 5)
plot_image_with_histogram(img_median_filtered, description="image égalisée et filtrée par médian")
Question 3 (1,5 points)¶
Pour la suite, vous allez avoir besoin d'utiliser le filtrage par convolution 2D. Implémentez la fonction conv2d qui prend en arguments une image et un masque et sort la convolution des deux. Vous n'utiliserez que les opérations matricielles et au maximum 2 boucles for, en particulier l'utilisation de toute fonction de convolution pré-existante est interdite. Pour que la taille de l'image convoluée soit la même que l'image de départ, vous effectuerez un 0 padding. Pour cela, vous pouvez vous inspirer du code suivant (à compléter). On supposera que le masque est toujours carré et sa taille toujours impaire.
def conv2d(img, mask):
out = np.zeros(img.shape, dtype=np.float32)
size_mask = mask.shape[0]
pad_values = ...
img = np.pad(img, ((pad_values, pad_values), (pad_values, pad_values)))
for i, row in enumerate(out):
for j, col in enumerate(row):
out[i, j] = ...
return out
Testez votre implémentation en filtrant votre image médiane avec le masque gaussien suivant:
mask = np.asarray([[1,2,1,2,1],[2,4,8,4,2],[1,8,18,8,1],[2,4,8,4,2],[1,2,1,2,1]])/90
et affichez le résultat. Qu'observez-vous?
def conv2d(img, mask):
assert mask.shape[0] == mask.shape[1], "mask must be a square matrix"
assert mask.shape[0] % 2 != 0, "kernel size must be odd"
out = np.zeros(img.shape, dtype=np.float32)
size_mask = mask.shape[0]
pad_values = size_mask // 2
img = np.pad(img, ((pad_values, pad_values), (pad_values, pad_values)))
flipped_mask = mask[::-1,::-1]
for i, row in enumerate(out):
for j, col in enumerate(row):
out[i, j] = (img[i: i+size_mask, j: j+size_mask] * flipped_mask).sum()
return out
mask = np.asarray([[1,2,1,2,1],[2,4,8,4,2],[1,8,18,8,1],[2,4,8,4,2],[1,2,1,2,1]])/90
convoluted_img = conv2d(img_median_filtered, mask)
plot_image_with_histogram(convoluted_img, description="image égalisée, filtrée et convoluée [après]")
plot_image_with_histogram(img_median_filtered, description="image égalisée et filtrée [avant]")
Réponse:
On remarque qu'en comparant l'image résultante avant et après l'application du filtre par la convolution que la nouvelle image est plus floue. On peut donc en déduire que le filtre était un filtre passe-bas, ce qui est confirmé en regardant le contenu du filtre, qui s'apparante bien à un filtre gaussien.
Question 4 (1 point)¶
Nous allons rectifier l'effet du filtre gaussien précédent en implémentant un filtre réhausseur de contours. Pour cela, nous ferons appel au Laplacien de l'image que l'on note: $\Delta(I)$. Implémentez un filtre réhausseur de contours fonctionnant en deux étapes:
- Un filtre passe-bas gaussien est appliqué: $I_g = G * I$
- $I_r = I_g + k\times \Delta(I_g)$
où $k$ est une constante que vous ferez varier pour obtenir le meilleur résultat.
On rappelle que le Laplacien peut s'obtenir avec une convolution avec un masque:Pour le filtre Gaussien, utilisez le filtre:
masque_laplacien = np.asarray([[-1,-1,-1], [-1,8,-1], [-1,-1,-1]])
masque_gaussien = np.asarray([[1,2,1], [2,4,2], [1,2,1]])
masque_laplacien = np.asarray([[-1,-1,-1], [-1,8,-1], [-1,-1,-1]])
masque_gaussien = np.asarray([[1,2,1], [2,4,2], [1,2,1]])/16
def enhance_edges(img, k, gaussian_mask, laplacien_mask):
img = img.astype(np.float32)
lp_img = cv2.filter2D(img, -1, kernel=gaussian_mask, borderType=cv2.BORDER_CONSTANT)
return lp_img + k * cv2.filter2D(lp_img, -1, kernel=laplacien_mask, borderType=cv2.BORDER_CONSTANT)
matplotlib.rcParams['figure.figsize'] = (15.0, 6.0)
for k in [0.1, 0.3, 0.5, 1, 3, 10]:
enhanced_edges_img = enhance_edges(img_median_filtered, k, masque_gaussien, masque_laplacien)
plot_image_with_histogram(enhanced_edges_img, description=f"image avec contours réhaussés et k={k}")
matplotlib.rcParams['figure.figsize'] = (15.0, 10.0)
Réponse: Le résultat qui semble avoir le contraste le plus élevé tout en évitant le bruit est celui où le paramètre k est fixé 0.3.
Exercice 2 - Compteur de pièces de monnaie¶
On vous a confié la tâche de réaliser un compteur de monnaie. Puisque la monnaie canadienne est essentiellement ronde, voilà l'occasion parfaite pour utiliser vos connaissances en filtrage morphologique.
Note: Il existe un excellent tutoriel sur l'utilisation du filtrage morphologique avec open-cv dont nous vous recommendons fortement la lecture.
Question 1 (0,5 point)¶
Chargez l'image pieces.jpg et convertissez la en niveau de gris. Affichez cette image.
FILENAME_EXO2 = 'pieces.jpg'
img = cv2.imread(FILENAME_EXO2, cv2.IMREAD_GRAYSCALE)
plt.imshow(img, vmin=0, vmax=255)
Question 2 (0,5 point)¶
Implémentez la fonction
def binariser(img, seuil):
...
qui met à 0 tous les pixels se trouvant en dessous du seuil, et à 255 tous les pixels se trouvant sur le seuil ou au dessus. Cette fonction prend en arguments une image et un seuil, puis retourne l'image binaire. Binarisez votre image avec un seuil égal à 250. Bien sûr, puisque nous voulons que ce soit les pièces qui soient blanches, vous devez inverser la couleur. Affichez l'image binaire.
def binariser(img, seuil, invert=True):
threshold_type = cv2.THRESH_BINARY_INV if invert else cv2.THRESH_BINARY
_, binary_img = cv2.threshold(img,seuil,255,threshold_type)
return binary_img.astype(np.bool_)
binary_img = binariser(img, 250, invert=True)
plt.imshow(binary_img, vmin=0, vmax=1)
Question 3 (1 point)¶
Vos pièces ne sont pas totalement pleines. Effectuez une fermeture sur l'image binaire afin de fermer les trous, tout en ne changeant pas trop la taille des pièces. Vous devez choisir un bon élément structurant. Affichez le résultat.
kernel = np.ones((11,11),np.uint8)
dilated_img = cv2.dilate(binary_img.astype(np.uint8),kernel,iterations = 1)
closing_img = cv2.erode(dilated_img,kernel,iterations = 1)
fig, ax = plt.subplots(ncols=2, figsize=(20, 10))
ax[0].imshow(binary_img, vmin=0, vmax=1)
ax[0].set_title("Avant la fermeture")
ax[1].imshow(closing_img, vmin=0, vmax=1)
ax[1].set_title("Après la fermeture");
fig.show();
Question 4 (3 points)¶
Implémentez la fonction:
def compter_monnaie(img_bin):
...
qui compte le total de la monnaie de votre image binaire à l'aide des opérateurs morphologiques imerode ou imdilate. La sortie de la fonction doit être exprimée en dollars. Voici des rayons de disques qui pourront vous être utiles : 200, 140, 120, 110, 90 (on ignorera les pièces de 1 dollar, 50 et 1 cent qui sont absentes de l'image).
Vous pouvez utiliser la fonction suivante pour compter le nombre de composantes connectées dans l'image.
def nb_components(img_bin):
num_labels, labels_im = cv2.connectedComponents(img_bin)
return num_labels-1 # On retire le fond qui est componente connectée mais qui ne nous intéresse pas
def nb_components(img_bin):
num_labels, labels_im = cv2.connectedComponents(img_bin)
return num_labels-1 # On retire le fond qui est componente connectée mais qui ne nous intéresse pas
def compter_monnaie(img_bin):
# ordered from biggest to smallest radius
coins_with_radius = [
(0.00, 200), # bottle cap
(2.00, 140),
(0.25, 120),
(0.05, 110),
(0.10, 90),
]
fig, ax = plt.subplots(ncols=len(coins_with_radius), figsize=(20, 40))
n_coins = np.zeros(len(coins_with_radius))
for idx, (coin, radius) in enumerate(coins_with_radius):
# find all cercles in which the kernel fits completely
kernel = cv2.circle(np.zeros((2*radius+1, 2*radius+1)), (radius, radius), radius, 1, -1).astype(np.uint8)
eroded_img = cv2.erode(closing_img.astype(np.uint8), kernel, iterations = 1)
ax[idx].imshow(eroded_img, vmin=0, vmax=1)
ax[idx].set_title(f'erode of coin {coin}$')
# substract the number of component of larger coins
n_coins[idx] = nb_components(eroded_img) - n_coins[:idx].sum()
plt.show()
print(*[f"There are {n_coin:.0f} coins of {coin_value:.2f}$." for n_coin, (coin_value, _) in zip(n_coins, coins_with_radius)], sep="\n")
return sum([n_coin * coin_value for n_coin, (coin_value, _) in zip(n_coins, coins_with_radius)])
monnaie = compter_monnaie(closing_img)
print(f"Total amount in the image: {monnaie:.2f}$")
Séance 2: Transformée de Fourier et Filtrage spectral¶
Exercice 3 Transformée de Fourier 2D¶
Cet exercice vise à vous familiariser avec le calcul de la transformée de Fourier en deux-dimensions et certaines de ses propriétés.
Pour cela, nous allons étudier la fonction: $f(x, y) = \sin(\frac{2\pi}{256}[f_1 x+f_2 y] ) \text{ avec } x, y \in [0, 256]$
Question 1 (1 point)¶
Créez la fonction:
def f(f1, f2):
...
qui prend en arguments les différentes constantes $\{f_1, f_2\}$ et renvoie une image monochrome de dimensions 256x256. Nous vous suggérons très fortement de vous référer à la documention de la fonction np.meshgrid pour une implémentation efficace (sans boucle for).
def f(f1, f2):
x = np.arange(256)
y = np.arange(256)
xx, yy = np.meshgrid(x, y)
return np.sin(2*np.pi*(1/256)*(f1 * xx + f2 * yy)).astype(np.uint8)
Question 2 (1,5 points)¶
Pour chaque jeux de paramètres suivants, affichez sur deux plots adjacents l'image de $f(x,y)$ et son spectre (amplitude de sa transformée de Fourier).
$\{f_1, f_2\}=$
- {12, 0}
- {0, 12}
- {12, 12}
- {12, 32}
- {32, -32}
Nous vous imposons de créer une fonction:
def fft_spectre(img):
...
qui prend en argument une image et renvoie le spectre de sa transformée de Fourier normalisée de telle sorte que sa valeur maximale soit égale à 1 et sa valeur minimale à 0. Utilisez les fonctions np.fft.fft2 et np.fft.fftshift. Pour la normalisation, vous pouvez utiliser la relation: $f_{normalisée} = \frac{f-\min(f)}{\max(f)-\min(f)}$ Note: N'oubliez pas de donner un titre à vos graphiques permettant de bien les identifier!
EPSILON = 1e-7 # Pour éviter de prendre le logarithme de zéro
def fft_spectre(img, scale):
spectrum = np.fft.fftshift(abs(np.fft.fft2(img)))
if scale == 'dB':
spectrum = 20*np.log10(spectrum + EPSILON)
minimum = np.min(spectrum)
maximum = np.max(spectrum)
normalized_spectrum = (spectrum - minimum) / (maximum - minimum) if maximum != minimum else spectrum
return normalized_spectrum
def plot_image_with_fft2(img, description, scale='linear'):
img = img.astype(np.uint8)
fig, ax = plt.subplots(ncols=2, figsize=(20, 8))
ax[0].imshow(img, vmin=0, vmax=255)
ax[0].set_title(f'{description}')
ax[1].imshow(fft_spectre(img, scale), vmin=0, vmax=1)
ax[1].set_title(f"Spectre d'amplitude de {description}");
fig.show()
params = [
[12, 0],
[0, 12],
[12, 12],
[12, 32],
[32, -32]
]
for f1, f2 in params:
img = f(f1, f2)
plot_image_with_fft2(img, f'Image générée avec f1 = {f1} et f2 = {f2}')
Question 3 (1 point)¶
Complétez les phrases suivantes, en vous basant sur vos observations à la question précédente sur les propriétés de la transformée de Fourier.
prop = "L'homothétie (agrandissement ou réduction) d'une image selon un axe donné \
entraîne %s sur son spectre dans le domaine de Fourier."
ans = "une homothétie inverse (réduction ou agrandissement)"
print(prop%ans)
prop = "La rotation d'angle alpha d'une image \
entraîne %s sur son spectre dans le domaine de Fourier."
ans = "une rotation d'angle alpha"
print(prop%ans)
Question 4 (1 point)¶
Soit la fonction:
$f(x, y)=\sin(\frac{2\pi}{256}f_1r)$ où $r=\sqrt{x^2+y^2} \text{ avec } x, y \in [-128, 128]$
Créez la fonction:
def wave(f1):
...
qui prend en arguments la constante $f_1$ et renvoie une image monochrome de dimensions 256x256.
Comme précédemment, tracez sur une figure le couple image/spectre pour les valeurs de $f_1$ suivantes: $f_1=\{12,64,128,256\}$ Que remarquez-vous?
def wave(f1):
x = np.arange(-128, 128)
y = np.arange(-128, 128)
xx, yy = np.meshgrid(x, y)
return np.sin(2*np.pi*(1/256)* f1 * np.sqrt(xx ** 2 + yy ** 2)).astype(np.uint8)
params = [12, 64, 128, 256]
for f in params:
img = wave(f)
plot_image_with_fft2(img, f'Image générée avec f = {f}')
Ce que l'on remarque¶
On peut voir qu'à patir d'une fréquence spatiale 128 Hz, il n'y a plus d'image qui est tracé. cela résulte en un spectre nul. Ceci est normal, car la fréquence de Nyquist pour ces images pour chacun des axes est de 128 Hz, soit la moitié de la résolution spatiale de l'image pour chacun des axes.
Exercice 4 Filtrage Spectral et homomorphique¶
Question 1 (1 point)¶
Un fichier mystère vous est fourni, représentant une célèbre image. Malheureusement, son identification est rendue impossible, du fait d'un double signal sinusoidal parasite, de fréquence $\pm 50$ selon chaque direction ($\pm x, \pm y$).
Chargez le fichier intitulé imageMystere.png, affichez-là ainsi que sa transformée de Fourier.
Note: Jusqu'à présent, nous avions étudié des images artificielles avec un contenu fréquentiel volontairement léger. Le spectre de l'image mystère est beaucoup plus dense, par conséquent sa visualisation est rendue légèrement plus complexe (les nombreuses fréquences présentes ont chacune des amplitudes très petites). Pour mieux voir ce qui se passe, affichez le spectre en décibel:
epsilon = 1e-7 # Pour éviter de prendre le logarithme de zéro
plt.imshow(20*np.log10(spectre+epsilon))
plt.show()
Sachant que le centre du spectre (après fftshift) correspond aux fréquences (0,0), vérifiez que vous retrouvez la trace des signaux pertubatoires sur le spectre de la fft.
FILENAME_EXO3 = 'imageMystere.png'
img = cv2.imread(FILENAME_EXO3, cv2.IMREAD_GRAYSCALE)
plot_image_with_fft2(img, 'Image mytérieuse', scale='dB')
Réponse: On voit effectivement quatres points blancs aux coordonnées (+-50, +-50) Hz. On peut également voir leurs harmoniques aux multiples de ces fréquences.
Question 2 (2 points)¶
Nous allons effacer les perturbations directement depuis le domaine spectrale, en créant un filtre qui sera directement appliqué (par multiplication) à la fft. Le filtre sera concu suivant un (ou plusieurs) profil(s) Gaussien selon la formule:
$H(u, v) = 1-e^{-((u \pm f)^2+(v \pm f)^2)/\sigma }$ où $\sigma$ permettra de contrôler la sélectivité du filtre et $f$ les fréquences à filtrer.
Adaptez cette formule pour créer un masque permettant de filtrer les perturbations et affichez ce masque. On choisira $\sigma=5$. Quel type de filtre (passe-bas, passe-haut, ...) reconnaissez-vous? Visuellement, que contrôle le paramètre $\sigma$ ?
def create_gaussian_profile(sigma, frequencies):
u = np.arange(-250, 250)
v = np.arange(-250, 250)
uu, vv = np.meshgrid(u, v)
h = np.ones((500, 500))
for freq_x, freq_y in frequencies:
h *= 1 - np.e ** (-((uu + freq_x)**2 + (vv + freq_y)**2)/sigma)
return h
# frequencies to filter : +- 50Hz in +- x and +- y
freqs = [
(-50, -50), (50, 50),
(50, -50), (-50, 50),
]
gaussian_profile = create_gaussian_profile(5, freqs)
fig, ax = plt.subplots(ncols=2, figsize=(20, 5))
ax[0].imshow(20*np.log10(gaussian_profile)+EPSILON)
ax[0].set_title('Filtre gaussien notch')
ax[1].imshow(fft_spectre(img, 'dB'), vmin=0, vmax=1)
ax[1].set_title("Spectre d'amplitude de l'image")
Type de filtre: On reconnait la signature spectrale de filtres spectraux locaux (notch filter). On voit deux paires de zones de rejet, en noir, qui sont symmétriques par rapport au centre du filtre.
Effet de sigma: De son nom et de son effet sur le filtre produit, nous pouvons voir que $\sigma$ représente l'écart type du gaussien. Ainsi, dans le cadre de la construction de filtre, un plus grand $\sigma$ signifie une plus grande zone de rejet, avec une pente moins abrupte, et inversement, un plus petit $\sigma$ signifie une petite zone de rejet avec une pente plus abrupte.
Question 3 (1 point)¶
Utilisez votre masque pour filtrer la transformée de Fourier (attention, pas juste le module!) et reconstruisez l'image par transformée de Fourier inverse np.fft.ifft2. Affichez l'image reconstruite ainsi que sa transformée de Fourier
Réponse: Nous remarquons encore des traces du bruit initial. Cela peut être dû par des harmoniques dans le spectre que nous n'avons pas filtré lors de l'étape précédente.
spectrum = np.fft.fft2(img)
spectrum_phase = np.angle(spectrum)
spectrum_amplitude = np.abs(spectrum)
spectrum_amplitude_centered = np.fft.fftshift(spectrum_amplitude)
spectrum_amplitude_centered_filtered = spectrum_amplitude_centered * gaussian_profile
filtered_img = np.fft.ifft2(
np.fft.ifftshift(spectrum_amplitude_centered_filtered) * np.exp(spectrum_phase*1j)
)
filtered_img = np.abs(filtered_img).astype(np.uint8)
plt.imshow(filtered_img)
Question 4 (1,5 points)¶
On rappelle la définition du filtrage homomorphique: $H(u, v)=(\gamma_H-\gamma_L)[1-e^{-c\frac{D^2(u, v)}{D^2_0}}]+\gamma_L$ où $D(u, v)=\sqrt{u^2+v^2}$
Implémentez ce filtre et testez-le avec les paramètres suivants: $D_0=2, \gamma_H=2, \gamma_L=0.5, c=1$. N'oubliez pas que ce filtre s'applique sur la transformée de Fourier de l'image logarithmique. Une fois le filtre appliqué et de retour dans le domaine spatial, il faudra donc prendre l'exponentielle du résultat.
Affichez le résultat (filtre et image filtrée) et décrivez l'effet des paramètres $\gamma_H, \gamma_L$ et $D_0$.
Réponse: Le paramètres ont les effets suivants: -$\gamma_H$ permet de définir le gain à appliquer aux hautes fréquences -$\gamma_L$ permet de définir l'atténuation à appliquer aux basses fréquences -$D_0$ permet de définir la distance à 0, ou la fréquence, à partir de laquelle appliquer le gain.
c = 1
def create_homomorphic_filter(gamma_h, gamma_l, d_0):
u = np.arange(-250, 250)
v = np.arange(-250, 250)
uu, vv = np.meshgrid(u, v)
distance = np.sqrt(uu ** 2 + vv ** 2)
exponent = -c * (distance**2 / (d_0 ** 2))
return (gamma_h - gamma_l)*(1 - np.e ** exponent) + gamma_l
homomorphic_filter = create_homomorphic_filter(**{'gamma_h': 2, 'gamma_l':0.5, 'd_0' : 2})
plt.imshow(homomorphic_filter)
log_filtered_img = np.log(filtered_img.astype(np.float64))
spectrum = np.fft.fft2(log_filtered_img)
spectrum_phase = np.angle(spectrum)
spectrum_amplitude = np.abs(spectrum)
spectrum_amplitude_centered = np.fft.fftshift(spectrum_amplitude)
spectrum_amplitude_centered_filtered = spectrum_amplitude_centered * homomorphic_filter
log_higlight_img = np.fft.ifft2(
np.fft.ifftshift(spectrum_amplitude_centered_filtered) * np.exp(spectrum_phase*1j)
)
higlight_img = np.exp(log_higlight_img)
plt.imshow(higlight_img.astype(np.uint8))
